Chapter 4 – Training Linear Models
This notebook contains all the sample code and solutions to the exercises in chapter 4.
本章では、最も単純なモデルのひとつである線形回帰モデルの詳細を見るところからスタートする。
大きく2つの異なる訓練方法を説明する。
1.閉形式の方程式...訓練セットに対してコスト関数が最小になるようなモデルパラメータを直接計算する。
2.反復的な最適化アプローチ...コスト関数が最小になるようにモデルパラメータを少しずつ操作し、収束させる勾配降下法。
Setup¶
First, let's make sure this notebook works well in both python 2 and 3, import a few common modules, ensure MatplotLib plots figures inline and prepare a function to save the figures:
# To support both python 2 and python 3
from __future__ import division, print_function, unicode_literals
# Common imports
import numpy as np
import os
# to make this notebook's output stable across runs
np.random.seed(42)
# To plot pretty figures
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)
# Where to save the figures
PROJECT_ROOT_DIR = "."
CHAPTER_ID = "training_linear_models"
def save_fig(fig_id, tight_layout=True):
path = os.path.join(PROJECT_ROOT_DIR, "images", CHAPTER_ID, fig_id + ".png")
print("Saving figure", fig_id)
if tight_layout:
plt.tight_layout()
plt.savefig(path, format='png', dpi=300)
# Ignore useless warnings (see SciPy issue #5998)
import warnings
warnings.filterwarnings(action="ignore", message="^internal gelsd")
4.1 線形回帰¶
回帰モデルの一般的な性能指標は、RMSE(第二章参照)であること、線形回帰モデルの訓練では、RMSEを最小にするパラメータΘを見つけることがゴールであることに留意する。
4.1.1 Linear regression using the Normal Equation (正規方程式)¶
コスト関数を最小にするΘの値を見つけるための閉形式解があり、言い換えれば直接結果を与えてくれるような数学的方程式であり、これを正規方程式という。
試しに、y=4+3x+ガウスノイズによって、テストデータを生成し、「正規方程式を使ってΘを計算し、モデルを生成 & プロット」してみる。
import numpy as np
import matplotlib.pyplot as plt
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
save_fig("generated_data_plot")
plt.show()
# 正規方程式
X_b = np.c_[np.ones((100, 1)), X] # add x0 = 1 to each instance
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
theta_best
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new] # add x0 = 1 to each instance
y_predict = X_new_b.dot(theta_best)
y_predict
plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.axis([0, 2, 0, 15])
plt.show()
The figure in the book actually corresponds to the following code, with a legend and axis labels:
plt.plot(X_new, y_predict, "r-", linewidth=2, label="Predictions")
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([0, 2, 0, 15])
save_fig("linear_model_predictions")
plt.show()
#numpyではなく、sklearnを使った場合のコード
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
lin_reg.intercept_, lin_reg.coef_
lin_reg.predict(X_new)
The LinearRegression class is based on the scipy.linalg.lstsq() function (the name stands for "least squares"), which you could call directly:
theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_b, y, rcond=1e-6)
theta_best_svd
This function computes $\mathbf{X}^+\mathbf{y}$, where $\mathbf{X}^{+}$ is the pseudoinverse of $\mathbf{X}$ (specifically the Moore-Penrose inverse). You can use np.linalg.pinv() to compute the pseudoinverse directly:
np.linalg.pinv(X_b).dot(y)
Note: the first releases of the book implied that the LinearRegression class was based on the Normal Equation. This was an error, my apologies: as explained above, it is based on the pseudoinverse, which ultimately relies on the SVD matrix decomposition of $\mathbf{X}$ (see chapter 8 for details about the SVD decomposition). Its time complexity is $O(n^2)$ and it works even when $m < n$ or when some features are linear combinations of other features (in these cases, $\mathbf{X}^T \mathbf{X}$ is not invertible so the Normal Equation fails), see issue #184 for more details. However, this does not change the rest of the description of the LinearRegression class, in particular, it is based on an analytical solution, it does not scale well with the number of features, it scales linearly with the number of instances, all the data must fit in memory, it does not require feature scaling and the order of the instances in the training set does not matter.
4.1.2 計算量¶
予想の対象インスタンスが2倍になったり、得著量が2倍になったりしても、予測にかかる時間はおおよそ2倍になる(=計算量は、予測したいインスタンス数、特徴量数に対して線形)
では、「特徴量が非常に多い」「訓練インスタンスがメモリに収まりきらない」ときに適している、別の線形回帰の訓練方法を見てみよう。
4.2 勾配降下法¶
基本的な考え方は、「コスト関数を最小にするために、パラメータを繰り返し操作する」。
操作...誤差関数の勾配を測定し、下降の方向にパラメータを調整する。
勾配降下法で重要な指標の一つに学習率というハイパーパラメータがある。この学習率が小さすぎると、収束までの反復数が増え、時間がかかる。逆に大きすぎると、別のパラメータ群を探索し始め、最悪発散して良いパラメータを見つけれられなくなる。
from IPython.display import Image
Image("./images/勾配降下法_1.PNG")
幸い、線形回帰モデルのMSEコスト関数は、凸関数(=全体の最小値が一つ)なので上記のように、「学習率を高めたため局所解に陥る」ということはない。
コラム¶
Image("./images/勾配降下法とスケーリング.jpg")
4.2.1 regression using batch gradient descent (バッチ勾配降下法)¶
勾配降下法を実装するには、個々のモデルパラメータΘjについて、コスト関数の勾配を計算する必要がある = Θjをほんのわずか変化させると、コスト関数がどれくらい変化するかを計算する必要がある(ここで偏微分が使われる)
コラム¶
個々のパラメータΘjをもとめる際、訓練セット全部に対する計算を行うアルゴリズムを特にバッチ勾配降下法という。
ステップごとに、訓練データ全体に対する勾配を計算するので、訓練セットが大規模だと遅くなる(※このデメリットを補った手法はこの後の節で説明する)
しかし、数十万個の特徴量がある線形回帰モデルの訓練では、正規方程式よりも、まだ勾配降下法を使った方がずっと高速である。
それでは、勾配降下法を実装してみよう
Image("./images/バッチ勾配降下法.PNG")
eta = 0.1 #学習率
n_iterations = 1000
m = 100
theta = np.random.randn(2,1) #ランダムな初期値
for iteration in range(n_iterations):
gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
theta = theta - eta * gradients
theta
結果は、正規方程式で解いた値と同じであることが見てとれる。
X_new_b.dot(theta)
別の学習率を与えた場合に、本手法を試した場合どうなるか。異なる3つの学習率を使った勾配降下法の最初の10ステップを見てみよう
theta_path_bgd = []
def plot_gradient_descent(theta, eta, theta_path=None):
m = len(X_b)
plt.plot(X, y, "b.")
n_iterations = 1000
for iteration in range(n_iterations):
if iteration < 10:
y_predict = X_new_b.dot(theta)
style = "b-" if iteration > 0 else "r--"
plt.plot(X_new, y_predict, style)
gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
theta = theta - eta * gradients
if theta_path is not None:
theta_path.append(theta)
plt.xlabel("$x_1$", fontsize=18)
plt.axis([0, 2, 0, 15])
plt.title(r"$\eta = {}$".format(eta), fontsize=16)
np.random.seed(42)
theta = np.random.randn(2,1) # random initialization
plt.figure(figsize=(10,4))
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)
save_fig("gradient_descent_plot")
plt.show()
左のグラフは学習率が低すぎて、最終的に解にたどり着くだろうが時間がかかる。
右は、学習率が高すぎて、アルゴリズムは発散し、求めたい解から離れている。
じゃあちょうどいい学習率を求めるにはどうすればいいかというと、2章で出てきたハイパーパラメータの調整の話(グリッドサーチとか)になる。
では例えばグリットサーチで求める際、なかなか収束しないモデルをグリットサーチでの探索候補から除く必要がある(めちゃくちゃ時間がかかるから)。 そのためには、反復回数を調整するというアプローチがある。この探索回数だが、小さすぎれば、まともに解を見つけられないままのモデルができ、逆に大きすぎれば、数回で解を実は見つけれれていて、パラメータの更新はもう起きないのに、無駄に勾配降下実行をする羽目になる。なので、シンプルなアプローチで、「反復回数はある程度大きくとるが、解の収束が見られ(=勾配ベクトルが小さくなる)たら、アルゴリズムを止めるという制御を加える」ことだ。
コラム¶
収束率...バッチ勾配降下法では、εの幅の中で最適解にたどり着くために、O(1/ε)回のイテレーションが必要となる = 許容誤差を1/10にすると、アルゴリズムの反復回数は10倍必要となる。
4.2.2 Stochastic Gradient Descent (確率的勾配降下法)¶
各ステップごとに、訓練セットからランダムに一つインスタンスを選び、それだけを使って、勾配を求める。一つのデータしか扱わないため、バッチ勾配降下法だとメモリに乗りきらないようなケース(データセットが強大)でも訓練できる。
ただ、ランダムにインスタンスを選ぶ性質上、コスト関数は、上下に動きならら(不規則に)収束していく。そのため、最終的に得られた結果は、バッチ勾配降下法で得られた解よりも悪いかもしれない。
Image("./images/SGD_2.PNG")
ただ、「ランダム性により局所解に陥る可能性は下がるが」「最小値に落ち着かないかもしれない」というジレンマがある。
対策としては、各イテレーションの学習率を決める関数(学習スケジュール)に従い、学習率を少しずつ下げて、全体の最小値を最終的に求められるようにする「焼きなまし法」というものがある。
以下に、簡単な学習スケジュールによる確率的勾配降下法の実装を示す。
theta_path_sgd = []
m = len(X_b) #m回のイテレーションを「1ラウンド」とカウントし、、各ラウンドを「エポック」と呼ぶ。
np.random.seed(42)
n_epochs = 50
t0, t1 = 5, 50 # learning schedule hyperparameters
def learning_schedule(t):
return t0 / (t + t1)
theta = np.random.randn(2,1) # random initialization
for epoch in range(n_epochs):
for i in range(m):
if epoch == 0 and i < 20: # not shown in the book
y_predict = X_new_b.dot(theta) # not shown
style = "b-" if i > 0 else "r--" # not shown
plt.plot(X_new, y_predict, style) # not shown
random_index = np.random.randint(m)
xi = X_b[random_index:random_index+1]
yi = y[random_index:random_index+1]
gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
eta = learning_schedule(epoch * m + i)
theta = theta - eta * gradients
theta_path_sgd.append(theta) # not shown
plt.plot(X, y, "b.") # not shown
plt.xlabel("$x_1$", fontsize=18) # not shown
plt.ylabel("$y$", rotation=0, fontsize=18) # not shown
plt.axis([0, 2, 0, 15]) # not shown
save_fig("sgd_plot") # not shown
plt.show() # not shown
theta
結果から、各ステップが不規則であるが、バッチ勾配降下法が1000回繰り替えしたのに対し、SGDは50回で同等の良い解にたどり着いている。
#sklearnで、SGDによる線形回帰を行う場合 (エポック= 50, 正規化なし、学習率 = 0.1)
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, tol=-np.infty, penalty=None, eta0=0.1, random_state=42)
sgd_reg.fit(X, y.ravel())
sgd_reg.intercept_, sgd_reg.coef_
4.2.3 Mini-batch gradient descent (ミニバッチ勾配降下法)¶
確率勾配降下法(SGD)が、各ステップで、1つのインスタンスを使うのに対し、複数のインスタンスをランダムに選んで使うのがミニバッチ勾配降下法。
ミニバッチがSGDよりも優れているのが、行列演算をGPUを使うことで速度向上が見込めるからだ。
theta_path_mgd = []
n_iterations = 50
minibatch_size = 20
np.random.seed(42)
theta = np.random.randn(2,1) # random initialization
t0, t1 = 200, 1000
def learning_schedule(t):
return t0 / (t + t1)
t = 0
for epoch in range(n_iterations):
shuffled_indices = np.random.permutation(m)
X_b_shuffled = X_b[shuffled_indices]
y_shuffled = y[shuffled_indices]
for i in range(0, m, minibatch_size):
t += 1
xi = X_b_shuffled[i:i+minibatch_size]
yi = y_shuffled[i:i+minibatch_size]
gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
eta = learning_schedule(t)
theta = theta - eta * gradients
theta_path_mgd.append(theta)
theta
theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)
plt.figure(figsize=(7,4))
plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], "r-s", linewidth=1, label="Stochastic")
plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], "g-+", linewidth=2, label="Mini-batch")
plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], "b-o", linewidth=3, label="Batch")
plt.legend(loc="upper left", fontsize=16)
plt.xlabel(r"$\theta_0$", fontsize=20)
plt.ylabel(r"$\theta_1$ ", fontsize=20, rotation=0)
plt.axis([2.5, 4.5, 2.3, 3.9])
save_fig("gradient_descent_paths_plot")
plt.show()
上記は、確率勾配降下法、ミニバッチ勾配降下法、バッチ勾配降下法それぞれのアルゴリズムにおける「パラメータ空間内の動き」を示している。
どれも最終的に最適解に近づくが、SGDはぴたっと最適解で止まるのに対し、他は変動が見られる。 ただし、SGDは実行時間が一番長く、確率勾配やミニバッチでも、適切な学習スケジュールを使えば最適解にたどり着けることを忘れてはならない。
4.3 Polynomial regression (多項式回帰)¶
データが単純な直線よりも難しい場合はどうだうか。
「非線形データに線形モデルを適合させること」ができ、シンプルなのは「多項式回帰」である。
多項式回帰...各特徴量の累乗を新しい特徴量として追加し、線形モデルで訓練する方法。
以下では、簡単な二次方程式から非線形データを生成し、多項式回帰でこれを予測してみる。
import numpy as np
import numpy.random as rnd
np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
save_fig("quadratic_data_plot")
plt.show()
明らかに、直線では予測できない。
では、各特徴量に二乗を新特徴量として訓練する。
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
X[0]
X_poly[0]
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
lin_reg.intercept_, lin_reg.coef_
X_new=np.linspace(-3, 3, 100).reshape(100, 1)
X_new_poly = poly_features.transform(X_new)
y_new = lin_reg.predict(X_new_poly)
plt.plot(X, y, "b.")
plt.plot(X_new, y_new, "r-", linewidth=2, label="Predictions")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([-3, 3, 0, 10])
save_fig("quadratic_predictions_plot")
plt.show()
なかなか良い予測である。
もとのデータは、y=0.5a^2 + 1.0b + 2.0に対して、モデルは y=0.56a^2 + 0.93b + 1.78という予測である。
複数の特徴量があるとき、多項式回帰は、「特徴量間の関係を見つけられる」というメリットがある。 なぜかというと、「PolynomialFeatures」関数は、指定された次数まで特徴量の全ての組み合わせを追加できるからである。 (ex.)aとbの特徴量があり、degree=2とした場合、a^2, b^2 だけでなく、組み合わせのabという特徴量も追加する)
4.4 学習曲線¶
高次の多項式を実行すれば、ただの線形回帰よりお訓練データにぴったりと適合させられる可能性が上がる。 下記に、さっきの訓練データに300次多項モデルを適用し、結果を純粋先鋭モデルと二次モデルの結果を比較する。
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
for style, width, degree in (("g-", 1, 300), ("b--", 2, 2), ("r-+", 2, 1)):
polybig_features = PolynomialFeatures(degree=degree, include_bias=False)
std_scaler = StandardScaler()
lin_reg = LinearRegression()
polynomial_regression = Pipeline([
("poly_features", polybig_features),
("std_scaler", std_scaler),
("lin_reg", lin_reg),
])
polynomial_regression.fit(X, y)
y_newbig = polynomial_regression.predict(X_new)
plt.plot(X_new, y_newbig, style, label=str(degree), linewidth=width)
plt.plot(X, y, "b.", linewidth=3)
plt.legend(loc="upper left")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
save_fig("high_degree_polynomials_plot")
plt.show()
もちろん、高次元多項式は、過学習しており、線形回帰モデルは過小適合となっている。 この場合、もっとも汎化するのはｍ二次モデルだが、「モデルをどの程度まで複雑にすべきか」「モデルがデータに過学習/過小適合していることをどうすれば見分けられるだろうか」。
2章では、「交差検証」により、モデルの汎化性能を推計した。モデルが訓練データに対して高い性能を発揮してもテストセットでうまく汎化していない場合は「過学習」であり、両方で性能が低ければ「過小適合」である。
もう一つの方法として、学習曲線の利用である。
学習曲線は、訓練セット(あるいは訓練イテレーション)、検証セットの性能をプロットしたもの。
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
def plot_learning_curves(model, X, y):
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=10)
train_errors, val_errors = [], []
for m in range(1, len(X_train)):
model.fit(X_train[:m], y_train[:m])
y_train_predict = model.predict(X_train[:m])
y_val_predict = model.predict(X_val)
train_errors.append(mean_squared_error(y_train[:m], y_train_predict))
val_errors.append(mean_squared_error(y_val, y_val_predict))
plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="train")
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="val")
plt.legend(loc="upper right", fontsize=14) # not shown in the book
plt.xlabel("Training set size", fontsize=14) # not shown
plt.ylabel("RMSE", fontsize=14) # not shown
# 線形回帰モデルの学習曲線
lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)
plt.axis([0, 80, 0, 3]) # not shown in the book
save_fig("underfitting_learning_curves_plot") # not shown
plt.show() # not shown
まず、訓練データに注目すると、訓練セットのインスタンスが1,2個なら、モデルは完全に適合するため、RMSEは0である。 ただ、訓練セットに新しいインスタンスが追加されると、ノイズやそもそも線形でないという理由から次第にモデルは訓練データに適合できなくなる。 そして、訓練誤差は次第に上昇し、やがてある地点で安定する。
次に検証セットｎ注目すると、少ないデータで訓練されたモデルでは十分に汎化できないため、最初は検証誤差が大きい(RMSE=2.4)が、データが増えると次第に適合し始める。しかし、直線ではデータを十分にモデリングできないため、検証誤差は下がらなくなり、ある地点で安定する。
これは過小適合モデルの典型的な例である。両方の曲線が一定の水準に達し、ともに非常に近接しているが、全体とて誤差RMSEが大きい。
コラム¶
モデルが訓練データに過小適合している場合は、訓練データを追加しても役に立たない。 より複雑なモデルを使うか、よりよい特徴量を用意する必要がある。
では、同じデータに対する10次元モデルの学習曲線を見てみよう。
from sklearn.pipeline import Pipeline
polynomial_regression = Pipeline([
("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
("lin_reg", LinearRegression()),
])
plot_learning_curves(polynomial_regression, X, y)
plt.axis([0, 80, 0, 3]) # not shown
save_fig("learning_curves_plot") # not shown
plt.show() # not shown
結果として、①だたの線形回帰より誤差が小さいが、②検証データに対する性能よりも訓練データに対する性能が「かなり高い」ため、過学習の特徴が見られる。しかし、訓練セットが上昇すると、両者の差は小さくなっているため、過学習度合も下がっているといえる。
コラム(バイアスと分散のトレードオフ)¶
統計学と機械学習の理論的研究から、モデルの汎化誤差は、3つの誤差の和として表現できる。
バイアス...でーたが実際には2次元なのに線形だと考えるなど、前提条件設定時の誤解である。(バイアスの高いモデルは過小適合しやすい)
分散...モデルが訓練データの小さな差異に過敏に反応するケースがある。自由度が高いモデルは分散が高くなりがちであり、訓練データに過学習する。
削減不能誤差...データ自体のノイズ。この部分の誤差を減らす唯一の方法は、データのクリーンアップなどである。
4.5 Regularized models(正則化された線形モデル)¶
モデルの正則化は、モデルの過学習を緩和するための有力な方法である。
線形回帰の正則化は、一般に、モデルの重みを制限して表現され、ここでは3種類の正則化方法を取り挙げる。
4.5.1 リッジ回帰¶
コスト関数に、正規化項として、αΣΘi^2を足して評価する。
これにより、学習アルゴリズムは、モデルの重さをできるだけ小さく保とうと学習する。 (正規化項は、訓練中のコスト関数にだけ加えられる)
このαが0のときは、ただの線形回帰になり、非常に大きい場合は全ての重みが限りなく0に近づく。
このリッジ回帰では、入力特徴量のスケールによって影響を受けるので、事前に、データスケーリングが必要となる。(ほとんどの正則化モデルにこれは当てはまる)
コラム¶
訓練中のコスト関数がテスト時に使う性能測定方法と異なることはよくある。
両者が異なるのは、「正規化」以外にも、「最適化しやすくするために導関数を持つよう訓練用コスト関数を用意」したり(例えばLog Lossなどのコスト関数で訓練し、適合率/再現率で評価するケース)
from sklearn.linear_model import Ridge
np.random.seed(42)
m = 20
X = 3 * np.random.rand(m, 1)
y = 1 + 0.5 * X + np.random.randn(m, 1) / 1.5
X_new = np.linspace(0, 3, 100).reshape(100, 1)
def plot_model(model_class, polynomial, alphas, **model_kargs):
for alpha, style in zip(alphas, ("b-", "g--", "r:")):
model = model_class(alpha, **model_kargs) if alpha > 0 else LinearRegression()
if polynomial:
model = Pipeline([
("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
("std_scaler", StandardScaler()),
("regul_reg", model),
])
model.fit(X, y)
y_new_regul = model.predict(X_new)
lw = 2 if alpha > 0 else 1
plt.plot(X_new, y_new_regul, style, linewidth=lw, label=r"$\alpha = {}$".format(alpha))
plt.plot(X, y, "b.", linewidth=3)
plt.legend(loc="upper left", fontsize=15)
plt.xlabel("$x_1$", fontsize=18)
plt.axis([0, 3, 0, 4])
plt.figure(figsize=(8,4))
plt.subplot(121)
plot_model(Ridge, polynomial=False, alphas=(0, 10, 100), random_state=42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Ridge, polynomial=True, alphas=(0, 10**-5, 1), random_state=42)
save_fig("ridge_regression_plot")
plt.show()
左は、プレーンなリッジモデルであり、線形予測である。
右は、「多項式回帰のときの要領でデータを拡張」「データスケーリング」した特徴量に対し、リッジモデルを適用したもの(リッジ正則化を伴った多項式回帰モデル)。
線形回帰の解き方として、閉形式、勾配降下法による解法があったが、リッジモデルも同様で、閉形式で計算してもいいし、勾配降下法で訓練してもよい。そのときの長所も短所も同じである。
#リッジモデルを閉形式で
from sklearn.linear_model import Ridge
ridge_reg = Ridge(alpha=1, solver="cholesky", random_state=42)
ridge_reg.fit(X, y)
ridge_reg.predict([[1.5]])
#リッジモデルを勾配降下法で
sgd_reg = SGDRegressor(max_iter=50, tol=-np.infty, penalty="l2", random_state=42)
sgd_reg.fit(X, y.ravel())
sgd_reg.predict([[1.5]])
ridge_reg = Ridge(alpha=1, solver="sag", random_state=42)
ridge_reg.fit(X, y)
ridge_reg.predict([[1.5]])
4.5.2 Lasso回帰¶
リッジ回帰と同様に正規化項を足すが、重みベクトルのL1ノルム(αΣ|Θi| )を使う。
from sklearn.linear_model import Lasso
plt.figure(figsize=(8,4))
plt.subplot(121)
plot_model(Lasso, polynomial=False, alphas=(0, 0.1, 1), random_state=42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Lasso, polynomial=True, alphas=(0, 10**-7, 1), tol=1, random_state=42)
save_fig("lasso_regression_plot")
plt.show()
Lasso回帰には、重要度の低い特徴量の重みを取り除く(つまり重みを0にする)傾向があるという重要な特徴がある。 例えば、上記右図の破線(α=1e - 07)はほとんど2次曲線に見える(=高次元多項式の重みの全てが0となっている。
言い換えれば、自動的に特徴量選択をし、疎なモデルを出力する。
Lassoコスト巻子は、Θi = 0で微分可能ではないが、Θi=0のき、代わりに劣勾配ベクトルを使えば、勾配降下法を適用できる。
次に示すのは、sklearnのLassoクラスを使った例である。
from sklearn.linear_model import Lasso
lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X, y)
lasso_reg.predict([[1.5]])
4.5.3 Elastic Net¶
正規化項は、リッジとLassoの混合であり、混ぜ方はハイパーパラメータrで任意設定できる。
from sklearn.linear_model import ElasticNet
elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5, random_state=42)
elastic_net.fit(X, y)
elastic_net.predict([[1.5]])
では、リッジ、Lasso、Elasticをどうつ使い分けるか。
Laasoは余計な特徴量の重みを0にしてくれるので、有力であるが、「訓練インタンス数 < 特徴量の数」、または「特徴量間に強い相関があるとき」、Lassoは不規則な動きを示す場合があるので、一般にLassoよりElastci Netの方がよい。
np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 2 + X + 0.5 * X**2 + np.random.randn(m, 1)
X_train, X_val, y_train, y_val = train_test_split(X[:50], y[:50].ravel(), test_size=0.5, random_state=10)
poly_scaler = Pipeline([
("poly_features", PolynomialFeatures(degree=90, include_bias=False)),
("std_scaler", StandardScaler()),
])
X_train_poly_scaled = poly_scaler.fit_transform(X_train)
X_val_poly_scaled = poly_scaler.transform(X_val)
sgd_reg = SGDRegressor(max_iter=1,
tol=-np.infty,
penalty=None,
eta0=0.0005,
warm_start=True,
learning_rate="constant",
random_state=42)
n_epochs = 500
train_errors, val_errors = [], []
for epoch in range(n_epochs):
sgd_reg.fit(X_train_poly_scaled, y_train)
y_train_predict = sgd_reg.predict(X_train_poly_scaled)
y_val_predict = sgd_reg.predict(X_val_poly_scaled)
train_errors.append(mean_squared_error(y_train, y_train_predict))
val_errors.append(mean_squared_error(y_val, y_val_predict))
best_epoch = np.argmin(val_errors)
best_val_rmse = np.sqrt(val_errors[best_epoch])
plt.annotate('Best model',
xy=(best_epoch, best_val_rmse),
xytext=(best_epoch, best_val_rmse + 1),
ha="center",
arrowprops=dict(facecolor='black', shrink=0.05),
fontsize=16,
)
best_val_rmse -= 0.03 # just to make the graph look better
plt.plot([0, n_epochs], [best_val_rmse, best_val_rmse], "k:", linewidth=2)
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="Validation set")
plt.plot(np.sqrt(train_errors), "r--", linewidth=2, label="Training set")
plt.legend(loc="upper right", fontsize=14)
plt.xlabel("Epoch", fontsize=14)
plt.ylabel("RMSE", fontsize=14)
save_fig("early_stopping_plot")
plt.show()
epoch(学習回数)が250を超えたあたりから、検証誤差が上がっている(=過学習し始める)。
早期打ち切りテクニックは、検証誤差が最小になったところで、訓練を中止する。
下記の早期打ち切りの基本実装は、「warm_start=True」となっているので、検証誤差が増加し始めたら、検証誤差が最小になったときにモデルパラメータにロールバックしてくれる。
from sklearn.base import clone
sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True, penalty=None,
learning_rate="constant", eta0=0.0005, random_state=42)
minimum_val_error = float("inf")
best_epoch = None
best_model = None
for epoch in range(1000):
sgd_reg.fit(X_train_poly_scaled, y_train) # continues where it left off
y_val_predict = sgd_reg.predict(X_val_poly_scaled)
val_error = mean_squared_error(y_val, y_val_predict)
if val_error < minimum_val_error:
minimum_val_error = val_error
best_epoch = epoch
best_model = clone(sgd_reg)
best_epoch, best_model
補足(Lasso正規化とリッジ正規化)¶
t1a, t1b, t2a, t2b = -1, 3, -1.5, 1.5
# ignoring bias term
t1s = np.linspace(t1a, t1b, 500)
t2s = np.linspace(t2a, t2b, 500)
t1, t2 = np.meshgrid(t1s, t2s)
T = np.c_[t1.ravel(), t2.ravel()]
Xr = np.array([[-1, 1], [-0.3, -1], [1, 0.1]])
yr = 2 * Xr[:, :1] + 0.5 * Xr[:, 1:]
J = (1/len(Xr) * np.sum((T.dot(Xr.T) - yr.T)**2, axis=1)).reshape(t1.shape)
N1 = np.linalg.norm(T, ord=1, axis=1).reshape(t1.shape)
N2 = np.linalg.norm(T, ord=2, axis=1).reshape(t1.shape)
t_min_idx = np.unravel_index(np.argmin(J), J.shape)
t1_min, t2_min = t1[t_min_idx], t2[t_min_idx]
t_init = np.array([[0.25], [-1]])
def bgd_path(theta, X, y, l1, l2, core = 1, eta = 0.1, n_iterations = 50):
path = [theta]
for iteration in range(n_iterations):
gradients = core * 2/len(X) * X.T.dot(X.dot(theta) - y) + l1 * np.sign(theta) + 2 * l2 * theta
theta = theta - eta * gradients
path.append(theta)
return np.array(path)
plt.figure(figsize=(12, 8))
for i, N, l1, l2, title in ((0, N1, 0.5, 0, "Lasso"), (1, N2, 0, 0.1, "Ridge")):
JR = J + l1 * N1 + l2 * N2**2
tr_min_idx = np.unravel_index(np.argmin(JR), JR.shape)
t1r_min, t2r_min = t1[tr_min_idx], t2[tr_min_idx]
levelsJ=(np.exp(np.linspace(0, 1, 20)) - 1) * (np.max(J) - np.min(J)) + np.min(J)
levelsJR=(np.exp(np.linspace(0, 1, 20)) - 1) * (np.max(JR) - np.min(JR)) + np.min(JR)
levelsN=np.linspace(0, np.max(N), 10)
path_J = bgd_path(t_init, Xr, yr, l1=0, l2=0)
path_JR = bgd_path(t_init, Xr, yr, l1, l2)
path_N = bgd_path(t_init, Xr, yr, np.sign(l1)/3, np.sign(l2), core=0)
plt.subplot(221 + i * 2)
plt.grid(True)
plt.axhline(y=0, color='k')
plt.axvline(x=0, color='k')
plt.contourf(t1, t2, J, levels=levelsJ, alpha=0.9)
plt.contour(t1, t2, N, levels=levelsN)
plt.plot(path_J[:, 0], path_J[:, 1], "w-o")
plt.plot(path_N[:, 0], path_N[:, 1], "y-^")
plt.plot(t1_min, t2_min, "rs")
plt.title(r"$\ell_{}$ penalty".format(i + 1), fontsize=16)
plt.axis([t1a, t1b, t2a, t2b])
if i == 1:
plt.xlabel(r"$\theta_1$", fontsize=20)
plt.ylabel(r"$\theta_2$", fontsize=20, rotation=0)
plt.subplot(222 + i * 2)
plt.grid(True)
plt.axhline(y=0, color='k')
plt.axvline(x=0, color='k')
plt.contourf(t1, t2, JR, levels=levelsJR, alpha=0.9)
plt.plot(path_JR[:, 0], path_JR[:, 1], "w-o")
plt.plot(t1r_min, t2r_min, "rs")
plt.title(title, fontsize=16)
plt.axis([t1a, t1b, t2a, t2b])
if i == 1:
plt.xlabel(r"$\theta_1$", fontsize=20)
save_fig("lasso_vs_ridge_plot")
plt.show()
4.6 Logistic regression(ロジスティック回帰)¶
回帰アルゴリズムのなかには、分類に使えるものがあり、ロジスティック回帰は、インスタンスが特定のクラスに属する確率を推計するためによく使われる。
4.6.1 確率の推計¶
ロジスティック回帰では、これまでの線形回帰モデルのように計算結果を直接するのではなく、0～1の間の値を出力する。
ロジスティック回帰モデルは、シグモイド関数(0～1の間の値を出力する)が利用される。そして実際の予測(2値分類)だと、出力結果が0.5より小さいならば、0, 大きければ1という予測により、クラス分類できる。
#シグモイド関数
t = np.linspace(-10, 10, 100)
sig = 1 / (1 + np.exp(-t))
plt.figure(figsize=(9, 3))
plt.plot([-10, 10], [0, 0], "k-")
plt.plot([-10, 10], [0.5, 0.5], "k:")
plt.plot([-10, 10], [1, 1], "k:")
plt.plot([0, 0], [-1.1, 1.1], "k-")
plt.plot(t, sig, "b-", linewidth=2, label=r"$\sigma(t) = \frac{1}{1 + e^{-t}}$")
plt.xlabel("t")
plt.legend(loc="upper left", fontsize=20)
plt.axis([-10, 10, -0.1, 1.1])
save_fig("logistic_function_plot")
plt.show()
4.6.2 訓練とコスト関数 (P.136)¶
どのように訓練するかだが、単一インスタンス用、訓練セット全体用にそれぞれコスト関数を分け、とりわけ訓練セット全体(全訓練インスタンスそれぞののコストの平均)は「Log Loss」と呼ばれる。
★ぶっちゃけなんで分けるか、理由がまだ調べきれていなので、あとで調べて展開します。
Log Lossを最小にするΘを計算する閉形式の方程式は存在しないが、勾配降下法で最小値を求めることはできる。
4.6.3 決定境界¶
2値分類において、陽インスタンスである確率と陰インスタンスがちょうど50%になる点。
例題として、3種類のiris(あやめ)を分類してみる。(データセットには、幅と長さの特徴量がある)
from sklearn import datasets
iris = datasets.load_iris()
list(iris.keys())
print(iris.DESCR)
幅の特徴量だけを用いて、ロジスティック回帰モデルを訓練&予測する。
X = iris["data"][:, 3:] # petal width
y = (iris["target"] == 2).astype(np.int) # 1 if Iris-Virginica, else 0
from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression(solver="liblinear", random_state=42)
log_reg.fit(X, y)
X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
plt.plot(X_new, y_proba[:, 1], "g-", linewidth=2, label="Iris-Virginica")
plt.plot(X_new, y_proba[:, 0], "b--", linewidth=2, label="Not Iris-Virginica")
The figure in the book actually is actually a bit fancier:
X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
decision_boundary = X_new[y_proba[:, 1] >= 0.5][0]
plt.figure(figsize=(8, 3))
plt.plot(X[y==0], y[y==0], "bs")
plt.plot(X[y==1], y[y==1], "g^")
plt.plot([decision_boundary, decision_boundary], [-1, 2], "k:", linewidth=2)
plt.plot(X_new, y_proba[:, 1], "g-", linewidth=2, label="Iris-Virginica")
plt.plot(X_new, y_proba[:, 0], "b--", linewidth=2, label="Not Iris-Virginica")
plt.text(decision_boundary+0.02, 0.15, "Decision boundary", fontsize=14, color="k", ha="center")
plt.arrow(decision_boundary, 0.08, -0.3, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')
plt.arrow(decision_boundary, 0.92, 0.3, 0, head_width=0.05, head_length=0.1, fc='g', ec='g')
plt.xlabel("Petal width (cm)", fontsize=14)
plt.ylabel("Probability", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 3, -0.02, 1.02])
save_fig("logistic_regression_plot")
plt.show()
decision_boundary
log_reg.predict([[1.7], [1.5]])
ちょうど幅(petal width)が1.6cmが決定協会となり、この値付近は、1.6cmよりもちょっとでも長ければバージニカ種、短かければ非バージニカ種と分類する。(=1.6cm付近の分類結果はあまり確証がもてない)
逆に、1cm未満だと、かなりの自信をもって(高い確率で)バージニカ種ではないと判断できる。
次に、幅だけでなく、長さの特徴量も加えて、同様に訓練&予測してみる。
from sklearn.linear_model import LogisticRegression
X = iris["data"][:, (2, 3)] # petal length, petal width
y = (iris["target"] == 2).astype(np.int)
log_reg = LogisticRegression(solver="liblinear", C=10**10, random_state=42)
log_reg.fit(X, y)
x0, x1 = np.meshgrid(
np.linspace(2.9, 7, 500).reshape(-1, 1),
np.linspace(0.8, 2.7, 200).reshape(-1, 1),
)
X_new = np.c_[x0.ravel(), x1.ravel()]
y_proba = log_reg.predict_proba(X_new)
plt.figure(figsize=(10, 4))
plt.plot(X[y==0, 0], X[y==0, 1], "bs")
plt.plot(X[y==1, 0], X[y==1, 1], "g^")
zz = y_proba[:, 1].reshape(x0.shape)
contour = plt.contour(x0, x1, zz, cmap=plt.cm.brg)
left_right = np.array([2.9, 7])
boundary = -(log_reg.coef_[0][0] * left_right + log_reg.intercept_[0]) / log_reg.coef_[0][1]
plt.clabel(contour, inline=1, fontsize=12)
plt.plot(left_right, boundary, "k--", linewidth=3)
plt.text(3.5, 1.5, "Not Iris-Virginica", fontsize=14, color="b", ha="center")
plt.text(6.5, 2.3, "Iris-Virginica", fontsize=14, color="g", ha="center")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.axis([2.9, 7, 0.8, 2.7])
save_fig("logistic_regression_contour_plot")
plt.show()
上記結果中の直線が、決定境界であり、モデルによれば右上の直線の向こう側の花のデータセットは、90%の確率でバージニカ種であり、左下の向こう側の花は、15%の確率でバージニカ種である、と分類している。
(澤田メモ)ロジスティック回帰において、特徴量が少ないと、50%の決定境界付近のデータ分類の精度がよくないから、特徴量を少なくしてはならないってことを言いたのかな
ソフトマックス回帰¶
ロジスティック回帰モデルは、3章の説明のゆおうに、複数の2項分類器を訓練して組み合わせなくても、複数のクラスを直接サポートできる。
それはソフトマックス回帰あるいは多項ロジスティック回帰と呼ぶ。
アルゴリズムは、ざっくり言うと、「個々のインスタンスにおいて、各クラスのスコア(下記画像何のUk)を計算して、ソフトマックス関数を適用」して、インスタンスがクラスkに属する確率を推計する。
コスト関数は交差エントロピー、出力はロジスティック回帰と同様に、推計された確率が最も高いクラスを予測として返す。
Image("./images/ソフトマックス回帰.PNG")
irisのあやめ3種類をクラス分類してみる。
X = iris["data"][:, (2, 3)] # petal length, petal width
y = iris["target"]
softmax_reg = LogisticRegression(multi_class="multinomial",solver="lbfgs", C=10, random_state=42)
softmax_reg.fit(X, y)
x0, x1 = np.meshgrid(
np.linspace(0, 8, 500).reshape(-1, 1),
np.linspace(0, 3.5, 200).reshape(-1, 1),
)
X_new = np.c_[x0.ravel(), x1.ravel()]
y_proba = softmax_reg.predict_proba(X_new)
y_predict = softmax_reg.predict(X_new)
zz1 = y_proba[:, 1].reshape(x0.shape)
zz = y_predict.reshape(x0.shape)
plt.figure(figsize=(10, 4))
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris-Virginica")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris-Versicolor")
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris-Setosa")
from matplotlib.colors import ListedColormap
custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])
plt.contourf(x0, x1, zz, cmap=custom_cmap)
contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)
plt.clabel(contour, inline=1, fontsize=12)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 7, 0, 3.5])
save_fig("softmax_regression_contour_plot")
plt.show()
softmax_reg.predict([[5, 2]])
softmax_reg.predict_proba([[5, 2]])
長さ5cm,幅2cmのあやめを見つけたときは、94%の確率でバージニカ種(クラス2)と予測した(あるいは、5.8%の確率でばバーシクル)。
Exercise solutions¶
1. to 11.¶
See appendix A.
12. Batch Gradient Descent with early stopping for Softmax Regression¶
(without using Scikit-Learn)
Let's start by loading the data. We will just reuse the Iris dataset we loaded earlier.
X = iris["data"][:, (2, 3)] # petal length, petal width
y = iris["target"]
We need to add the bias term for every instance ($x_0 = 1$):
X_with_bias = np.c_[np.ones([len(X), 1]), X]
And let's set the random seed so the output of this exercise solution is reproducible:
np.random.seed(2042)
The easiest option to split the dataset into a training set, a validation set and a test set would be to use Scikit-Learn's train_test_split() function, but the point of this exercise is to try understand the algorithms by implementing them manually. So here is one possible implementation:
test_ratio = 0.2
validation_ratio = 0.2
total_size = len(X_with_bias)
test_size = int(total_size * test_ratio)
validation_size = int(total_size * validation_ratio)
train_size = total_size - test_size - validation_size
rnd_indices = np.random.permutation(total_size)
X_train = X_with_bias[rnd_indices[:train_size]]
y_train = y[rnd_indices[:train_size]]
X_valid = X_with_bias[rnd_indices[train_size:-test_size]]
y_valid = y[rnd_indices[train_size:-test_size]]
X_test = X_with_bias[rnd_indices[-test_size:]]
y_test = y[rnd_indices[-test_size:]]
The targets are currently class indices (0, 1 or 2), but we need target class probabilities to train the Softmax Regression model. Each instance will have target class probabilities equal to 0.0 for all classes except for the target class which will have a probability of 1.0 (in other words, the vector of class probabilities for ay given instance is a one-hot vector). Let's write a small function to convert the vector of class indices into a matrix containing a one-hot vector for each instance:
def to_one_hot(y):
n_classes = y.max() + 1
m = len(y)
Y_one_hot = np.zeros((m, n_classes))
Y_one_hot[np.arange(m), y] = 1
return Y_one_hot
Let's test this function on the first 10 instances:
y_train[:10]
to_one_hot(y_train[:10])
Looks good, so let's create the target class probabilities matrix for the training set and the test set:
Y_train_one_hot = to_one_hot(y_train)
Y_valid_one_hot = to_one_hot(y_valid)
Y_test_one_hot = to_one_hot(y_test)
Now let's implement the Softmax function. Recall that it is defined by the following equation:
$\sigma\left(\mathbf{s}(\mathbf{x})\right)_k = \dfrac{\exp\left(s_k(\mathbf{x})\right)}{\sum\limits_{j=1}^{K}{\exp\left(s_j(\mathbf{x})\right)}}$
def softmax(logits):
exps = np.exp(logits)
exp_sums = np.sum(exps, axis=1, keepdims=True)
return exps / exp_sums
We are almost ready to start training. Let's define the number of inputs and outputs:
n_inputs = X_train.shape[1] # == 3 (2 features plus the bias term)
n_outputs = len(np.unique(y_train)) # == 3 (3 iris classes)
Now here comes the hardest part: training! Theoretically, it's simple: it's just a matter of translating the math equations into Python code. But in practice, it can be quite tricky: in particular, it's easy to mix up the order of the terms, or the indices. You can even end up with code that looks like it's working but is actually not computing exactly the right thing. When unsure, you should write down the shape of each term in the equation and make sure the corresponding terms in your code match closely. It can also help to evaluate each term independently and print them out. The good news it that you won't have to do this everyday, since all this is well implemented by Scikit-Learn, but it will help you understand what's going on under the hood.
So the equations we will need are the cost function:
$J(\mathbf{\Theta}) =
- \dfrac{1}{m}\sum\limits{i=1}^{m}\sum\limits{k=1}^{K}{y_k^{(i)}\log\left(\hat{p}_k^{(i)}\right)}$
And the equation for the gradients:
$\nabla_{\mathbf{\theta}^{(k)}} \, J(\mathbf{\Theta}) = \dfrac{1}{m} \sum\limits_{i=1}^{m}{ \left ( \hat{p}^{(i)}_k - y_k^{(i)} \right ) \mathbf{x}^{(i)}}$
Note that $\log\left(\hat{p}_k^{(i)}\right)$ may not be computable if $\hat{p}_k^{(i)} = 0$. So we will add a tiny value $\epsilon$ to $\log\left(\hat{p}_k^{(i)}\right)$ to avoid getting nan values.
eta = 0.01
n_iterations = 5001
m = len(X_train)
epsilon = 1e-7
Theta = np.random.randn(n_inputs, n_outputs)
for iteration in range(n_iterations):
logits = X_train.dot(Theta)
Y_proba = softmax(logits)
loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon), axis=1))
error = Y_proba - Y_train_one_hot
if iteration % 500 == 0:
print(iteration, loss)
gradients = 1/m * X_train.T.dot(error)
Theta = Theta - eta * gradients
And that's it! The Softmax model is trained. Let's look at the model parameters:
Theta
Let's make predictions for the validation set and check the accuracy score:
logits = X_valid.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba, axis=1)
accuracy_score = np.mean(y_predict == y_valid)
accuracy_score
Well, this model looks pretty good. For the sake of the exercise, let's add a bit of $\ell_2$ regularization. The following training code is similar to the one above, but the loss now has an additional $\ell_2$ penalty, and the gradients have the proper additional term (note that we don't regularize the first element of Theta since this corresponds to the bias term). Also, let's try increasing the learning rate eta.
eta = 0.1
n_iterations = 5001
m = len(X_train)
epsilon = 1e-7
alpha = 0.1 # regularization hyperparameter
Theta = np.random.randn(n_inputs, n_outputs)
for iteration in range(n_iterations):
logits = X_train.dot(Theta)
Y_proba = softmax(logits)
xentropy_loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon), axis=1))
l2_loss = 1/2 * np.sum(np.square(Theta[1:]))
loss = xentropy_loss + alpha * l2_loss
error = Y_proba - Y_train_one_hot
if iteration % 500 == 0:
print(iteration, loss)
gradients = 1/m * X_train.T.dot(error) + np.r_[np.zeros([1, n_outputs]), alpha * Theta[1:]]
Theta = Theta - eta * gradients
Because of the additional $\ell_2$ penalty, the loss seems greater than earlier, but perhaps this model will perform better? Let's find out:
logits = X_valid.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba, axis=1)
accuracy_score = np.mean(y_predict == y_valid)
accuracy_score
Cool, perfect accuracy! We probably just got lucky with this validation set, but still, it's pleasant.
Now let's add early stopping. For this we just need to measure the loss on the validation set at every iteration and stop when the error starts growing.
eta = 0.1
n_iterations = 5001
m = len(X_train)
epsilon = 1e-7
alpha = 0.1 # regularization hyperparameter
best_loss = np.infty
Theta = np.random.randn(n_inputs, n_outputs)
for iteration in range(n_iterations):
logits = X_train.dot(Theta)
Y_proba = softmax(logits)
xentropy_loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon), axis=1))
l2_loss = 1/2 * np.sum(np.square(Theta[1:]))
loss = xentropy_loss + alpha * l2_loss
error = Y_proba - Y_train_one_hot
gradients = 1/m * X_train.T.dot(error) + np.r_[np.zeros([1, n_outputs]), alpha * Theta[1:]]
Theta = Theta - eta * gradients
logits = X_valid.dot(Theta)
Y_proba = softmax(logits)
xentropy_loss = -np.mean(np.sum(Y_valid_one_hot * np.log(Y_proba + epsilon), axis=1))
l2_loss = 1/2 * np.sum(np.square(Theta[1:]))
loss = xentropy_loss + alpha * l2_loss
if iteration % 500 == 0:
print(iteration, loss)
if loss < best_loss:
best_loss = loss
else:
print(iteration - 1, best_loss)
print(iteration, loss, "early stopping!")
break
logits = X_valid.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba, axis=1)
accuracy_score = np.mean(y_predict == y_valid)
accuracy_score
Still perfect, but faster.
Now let's plot the model's predictions on the whole dataset:
x0, x1 = np.meshgrid(
np.linspace(0, 8, 500).reshape(-1, 1),
np.linspace(0, 3.5, 200).reshape(-1, 1),
)
X_new = np.c_[x0.ravel(), x1.ravel()]
X_new_with_bias = np.c_[np.ones([len(X_new), 1]), X_new]
logits = X_new_with_bias.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba, axis=1)
zz1 = Y_proba[:, 1].reshape(x0.shape)
zz = y_predict.reshape(x0.shape)
plt.figure(figsize=(10, 4))
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris-Virginica")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris-Versicolor")
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris-Setosa")
from matplotlib.colors import ListedColormap
custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])
plt.contourf(x0, x1, zz, cmap=custom_cmap)
contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)
plt.clabel(contour, inline=1, fontsize=12)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="upper left", fontsize=14)
plt.axis([0, 7, 0, 3.5])
plt.show()
And now let's measure the final model's accuracy on the test set:
logits = X_test.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba, axis=1)
accuracy_score = np.mean(y_predict == y_test)
accuracy_score
Our perfect model turns out to have slight imperfections. This variability is likely due to the very small size of the dataset: depending on how you sample the training set, validation set and the test set, you can get quite different results. Try changing the random seed and running the code again a few times, you will see that the results will vary.